import numpy as np
from scipy.sparse import diags, eye, bmat, csc_matrix
from scipy.sparse.linalg import spsolve

def solve_poisson_equation_direct(Nx, Ny):
    """
    Solves the Poisson equation using direct method.
    
    Args:
        Nx (int): Number of grid points in x-direction.
        Ny (int): Number of grid points in y-direction.
        
    Returns:
        tuple: (x, y, u), where x, y are the coordinate grids and u is the solution matrix.
    """
    assert isinstance(Nx, int) and Nx > 1, "Nx must be an integer greater than 1."
    assert isinstance(Ny, int) and Ny > 1, "Ny must be an integer greater than 1."
    
    # Grid spacing
    dx = 2 / Nx
    dy = 2 / Ny
    
    # Grid generation
    xx = np.linspace(-1 + dx/2, 1 - dx/2, Nx)
    yy = np.linspace(-1 + dy/2, 1 - dy/2, Ny)
    x, y = np.meshgrid(xx, yy)
    
    # Construction of A matrix for x derivative
    tmp = -2 * np.ones(Nx)
    tmp[[0, -1]] = -1  # Adjusting ends for Neumann BC
    Ax = diags([tmp[1:], 1, tmp[:-1]], [-1, 0, 1]).toarray()  # Convert to array for clarity
    
    # Block diagonal construction for y derivative
    dd = np.eye(Nx)
    ca = Ax / dx**2 - 2 * dd / dy**2
    sca = csc_matrix(ca)  # Sparse representation
    blocks = [sca for _ in range(Ny)]
    blocks[0] = csc_matrix(Ax / dx**2 - dd / dy**2)  # Adjust first block for Neumann BC
    blocks[-1] = blocks[0]  # Adjust last block similarly
    Abig = bmat(blocks, format='csc')
    
    # Modify Abig to include y-derivative interactions
    d4y = eye((Ny-1)*Nx, format='csc')
    Abig = Abig + csc_matrix((d4y / dy**2, (range(Nx), range(Nx, (Ny-1)*Nx))), shape=Abig.shape)
    Abig = Abig + csc_matrix((d4y / dy**2, (range(Nx, Abig.shape[0]), range(Nx))))

    # Make Abig singular matrix non-singular by pinning a solution point
    Abig[0, :] = 0
    Abig[0, 0] = 1
    b = np.ravel(-2 * np.pi**2 * np.cos(np.pi*x) * np.cos(np.pi*y))  # Right-hand side vector
    b[0] = 0
    
    # Solve the linear system
    u = spsolve(Abig, b)
    u = u.reshape(Nx, Ny)
    
    return x, y, u